Chapter 6 Beta diversity

load("data/data.Rdata")

beta_q0n <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  dplyr::select_if(~!all(. == 0)) %>%
  hillpair(., q = 0)

beta_q1n <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  dplyr::select_if(~!all(. == 0)) %>%
  hillpair(., q = 1)

beta_q1p <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  dplyr::select_if(~!all(. == 0)) %>%
  hillpair(., q = 1, tree = genome_tree)

beta_q1f <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  dplyr::select_if(~!all(. == 0)) %>%
  hillpair(., q = 1, dist = dist)

6.1 location

6.1.1 Richness diversity

betadisper(beta_q0n$S, sample_metadata$location) %>% permutest(., pairwise=TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)    
Groups     5 0.51463 0.102926 11.076    999  0.001 ***
Residuals 86 0.79917 0.009293                         
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
               Aruba     Brazil  CaboVerde      Spain    Denmark Malaysia
Aruba                1.0000e-03 2.0000e-03 1.0000e-03 1.0000e-03    0.015
Brazil    8.1559e-05            1.3700e-01 6.8000e-02 9.0600e-01    0.151
CaboVerde 6.2922e-04 1.1254e-01            1.0000e-03 3.0000e-02    0.937
Spain     1.7363e-12 7.3505e-02 3.5940e-05            2.1000e-02    0.001
Denmark   1.3206e-07 9.0556e-01 2.9720e-02 2.3526e-02               0.051
Malaysia  9.6128e-03 1.3730e-01 9.2726e-01 4.6570e-04 5.9273e-02         

adonis2(beta_q0n$S ~ location, 
        data = sample_metadata %>% arrange(match(sample,labels(beta_q0n$S))), 
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
location	5	6.57607	0.2442751	5.559605	0.001
Residual	86	20.34468	0.7557249	NA	NA
Total	91	26.92075	1.0000000	NA	NA

pairwise.adonis(beta_q0n$S, sample_metadata$location, perm = 999)

                   pairs Df SumsOfSqs  F.Model         R2 p.value p.adjusted sig
1        Aruba vs Brazil  1 2.0520609 7.280171 0.19528265   0.001      0.015   .
2     Aruba vs CaboVerde  1 1.0026264 3.213494 0.09975614   0.002      0.030   .
3       Aruba vs Denmark  1 1.9513181 7.018855 0.19486614   0.001      0.015   .
4      Aruba vs Malaysia  1 1.1650328 3.678697 0.11257172   0.001      0.015   .
5         Aruba vs Spain  1 2.0433431 8.167535 0.21974917   0.001      0.015   .
6    Brazil vs CaboVerde  1 2.1090004 9.164611 0.24013374   0.001      0.015   .
7      Brazil vs Denmark  1 0.5703370 2.907951 0.09113563   0.001      0.015   .
8     Brazil vs Malaysia  1 0.9283278 3.953419 0.11996993   0.001      0.015   .
9        Brazil vs Spain  1 0.6343957 3.769479 0.11503016   0.001      0.015   .
10  CaboVerde vs Denmark  1 1.8099512 8.070070 0.22373314   0.001      0.015   .
11 CaboVerde vs Malaysia  1 1.2143821 4.593887 0.14094322   0.001      0.015   .
12    CaboVerde vs Spain  1 1.8141430 9.281722 0.24896172   0.001      0.015   .
13   Denmark vs Malaysia  1 1.0124763 4.418610 0.13629855   0.001      0.015   .
14      Denmark vs Spain  1 0.5305346 3.310771 0.10573905   0.001      0.015   .
15     Malaysia vs Spain  1 0.8450695 4.218747 0.13094076   0.001      0.015   .

#pdf("figures/beta_q0_loca.pdf",width=9, height=5)
beta_q0n$S %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = "sample") %>%
  group_by(location) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = location, fill = location)) +
    geom_point(size = 4) +
    #   stat_ellipse(aes(color = beta_q1n_nmds$Groups))+
    geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
    scale_color_manual(values = location_colors)+
    theme_classic() +
    theme(
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      axis.title = element_text(size = 20, face = "bold"),
      axis.text = element_text(face = "bold", size = 18),
      panel.background = element_blank(),
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      legend.text = element_text(size = 16),
      legend.title = element_text(size = 18),
      legend.position = "right", legend.box = "vertical"
    ) 

#dev.off()

6.1.2 Neutral diversity

betadisper(beta_q1n$S, sample_metadata$location) %>% permutest(., pairwise=TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df Sum Sq  Mean Sq      F N.Perm Pr(>F)   
Groups     5 0.2728 0.054560 3.7079    999  0.006 **
Residuals 86 1.2654 0.014715                        
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
               Aruba     Brazil  CaboVerde      Spain    Denmark Malaysia
Aruba                1.1000e-02 1.0000e-03 1.0000e-03 1.0000e-03    0.002
Brazil    9.2383e-03            7.7300e-01 2.4400e-01 6.4900e-01    0.881
CaboVerde 7.7521e-04 7.7033e-01            3.1300e-01 8.8200e-01    0.890
Spain     1.6646e-05 2.2327e-01 3.0424e-01            3.7200e-01    0.258
Denmark   3.3254e-04 6.6322e-01 8.7642e-01 3.7291e-01               0.785
Malaysia  1.1422e-03 8.5989e-01 8.9763e-01 2.4754e-01 7.7435e-01         

adonis2(beta_q1n$S ~ location, 
        data = sample_metadata %>% arrange(match(sample,labels(beta_q1n$S))), 
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
location	5	5.709562	0.221304	4.88821	0.001
Residual	86	20.090068	0.778696	NA	NA
Total	91	25.799630	1.000000	NA	NA

pairwise.adonis(beta_q1n$S, sample_metadata$location, perm = 999)

                   pairs Df SumsOfSqs  F.Model         R2 p.value p.adjusted sig
1        Aruba vs Brazil  1 1.3739907 4.789154 0.13766227   0.001      0.015   .
2     Aruba vs CaboVerde  1 1.0636589 3.806728 0.11603497   0.001      0.015   .
3       Aruba vs Denmark  1 1.6037782 5.813780 0.16699651   0.001      0.015   .
4      Aruba vs Malaysia  1 1.2584980 4.466640 0.13346545   0.001      0.015   .
5         Aruba vs Spain  1 1.4382789 5.558030 0.16083181   0.001      0.015   .
6    Brazil vs CaboVerde  1 1.5901673 7.027154 0.19505161   0.001      0.015   .
7      Brazil vs Denmark  1 0.6954159 3.122217 0.09719804   0.001      0.015   .
8     Brazil vs Malaysia  1 0.7313843 3.199010 0.09935119   0.001      0.015   .
9        Brazil vs Spain  1 0.4007935 1.948927 0.06297236   0.034      0.510    
10  CaboVerde vs Denmark  1 1.8198218 8.556140 0.23405480   0.001      0.015   .
11 CaboVerde vs Malaysia  1 1.3171418 6.019861 0.17695136   0.001      0.015   .
12    CaboVerde vs Spain  1 1.5415454 7.905417 0.22017337   0.001      0.015   .
13   Denmark vs Malaysia  1 1.0691578 4.970172 0.15074753   0.001      0.015   .
14      Denmark vs Spain  1 0.5766622 3.014217 0.09718823   0.004      0.060    
15     Malaysia vs Spain  1 0.6446972 3.265586 0.10444666   0.001      0.015   .

#pdf("figures/beta_q1n_loca.pdf",width=9, height=5)
beta_q1n$S %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = "sample") %>%
  group_by(location) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = location, fill = location)) +
    geom_point(size = 4) +
    #   stat_ellipse(aes(color = beta_q1n_nmds$Groups))+
    geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
    scale_color_manual(values = location_colors)+
    theme_classic() +
    theme(
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      axis.title = element_text(size = 20, face = "bold"),
      axis.text = element_text(face = "bold", size = 18),
      panel.background = element_blank(),
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      legend.text = element_text(size = 16),
      legend.title = element_text(size = 18),
      legend.position = "right", legend.box = "vertical"
    ) 

#dev.off()

6.1.3 Phylogenetic diversity

betadisper(beta_q1p$S, sample_metadata$location) %>% permutest(., pairwise=TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df Sum Sq Mean Sq      F N.Perm Pr(>F)  
Groups     5 0.2152 0.04304 2.1047    999  0.082 .
Residuals 86 1.7587 0.02045                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
               Aruba     Brazil  CaboVerde      Spain    Denmark Malaysia
Aruba                0.06500000 0.04800000 0.00200000 0.10900000    0.004
Brazil    0.06236518            0.85700000 0.33200000 0.66800000    0.679
CaboVerde 0.05279004 0.86912737            0.20700000 0.74600000    0.515
Spain     0.00060645 0.32391452 0.20693940            0.11000000    0.480
Denmark   0.09811567 0.63811338 0.73976492 0.10108967               0.300
Malaysia  0.00519311 0.67845753 0.52055144 0.50696000 0.31212601         

adonis2(beta_q1p$S ~ location, 
        data = sample_metadata %>% arrange(match(sample,labels(beta_q1p$S))), 
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
location	5	2.690791	0.2683828	6.309563	0.001
Residual	86	7.335152	0.7316172	NA	NA
Total	91	10.025944	1.0000000	NA	NA

pairwise.adonis(beta_q1p$S, sample_metadata$location, perm = 999)

                   pairs Df  SumsOfSqs    F.Model         R2 p.value p.adjusted sig
1        Aruba vs Brazil  1 0.64714035  6.1133065 0.16928127   0.001      0.015   .
2     Aruba vs CaboVerde  1 0.56229778  5.3504131 0.15575979   0.001      0.015   .
3       Aruba vs Denmark  1 1.03758657  9.5591181 0.24790811   0.001      0.015   .
4      Aruba vs Malaysia  1 0.48047063  5.0352414 0.14794199   0.001      0.015   .
5         Aruba vs Spain  1 0.60047832  6.7919744 0.18976250   0.001      0.015   .
6    Brazil vs CaboVerde  1 0.70961443  8.1023431 0.21837821   0.001      0.015   .
7      Brazil vs Denmark  1 0.32742702  3.5968626 0.11034383   0.014      0.210    
8     Brazil vs Malaysia  1 0.19570506  2.5119802 0.07971509   0.031      0.465    
9        Brazil vs Spain  1 0.06456161  0.9106384 0.03044530   0.435      1.000    
10  CaboVerde vs Denmark  1 1.07120011 11.9405505 0.29895808   0.001      0.015   .
11 CaboVerde vs Malaysia  1 0.50032095  6.5728056 0.19011490   0.002      0.030   .
12    CaboVerde vs Spain  1 0.82104026 11.9236999 0.29866220   0.001      0.015   .
13   Denmark vs Malaysia  1 0.54258755  6.8084757 0.19559821   0.001      0.015   .
14      Denmark vs Spain  1 0.24626102  3.3999417 0.10827860   0.017      0.255    
15     Malaysia vs Spain  1 0.26932060  4.5771932 0.14050299   0.001      0.015   .

#pdf("figures/beta_q1p_loca.pdf",width=9, height=5)
beta_q1n$S %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = "sample") %>%
  group_by(location) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = location, fill = location)) +
    geom_point(size = 4) +
    #   stat_ellipse(aes(color = beta_q1n_nmds$Groups))+
    geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
    scale_color_manual(values = location_colors)+
    theme_classic() +
    theme(
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      axis.title = element_text(size = 20, face = "bold"),
      axis.text = element_text(face = "bold", size = 18),
      panel.background = element_blank(),
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      legend.text = element_text(size = 16),
      legend.title = element_text(size = 18),
      legend.position = "right", legend.box = "vertical"
    ) 

#dev.off()

6.1.4 Functional diversity

betadisper(beta_q1f$S, sample_metadata$location) %>% permutest(., pairwise=TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df Sum Sq  Mean Sq      F N.Perm Pr(>F)    
Groups     5 1.5016 0.300325 9.1732    999  0.001 ***
Residuals 86 2.8156 0.032739                         
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
               Aruba     Brazil  CaboVerde      Spain    Denmark Malaysia
Aruba                3.0000e-03 6.0600e-01 3.0000e-03 1.0000e-03    0.497
Brazil    8.6955e-04            2.0000e-03 6.9600e-01 2.4400e-01    0.003
CaboVerde 5.8936e-01 9.7727e-04            3.0000e-03 1.0000e-03    0.873
Spain     1.5125e-03 6.8358e-01 1.5804e-03            3.4000e-02    0.001
Denmark   1.7539e-04 2.3957e-01 9.2505e-05 3.3658e-02               0.001
Malaysia  4.9615e-01 1.0071e-03 8.8292e-01 1.5959e-03 7.4748e-05         

adonis2(beta_q1f$S ~ location, 
        data = sample_metadata %>% arrange(match(sample,labels(beta_q1f$S))), 
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
location	5	8.815599	0.5844961	24.19552	0.001
Residual	86	6.266793	0.4155039	NA	NA
Total	91	15.082391	1.0000000	NA	NA

pairwise.adonis(beta_q1f$S, sample_metadata$location, perm = 999)

                   pairs Df   SumsOfSqs    F.Model           R2 p.value p.adjusted sig
1        Aruba vs Brazil  1  3.77109355 40.6261637  0.575228238   0.001      0.015   .
2     Aruba vs CaboVerde  1  0.81599223  5.5801053  0.161367505   0.013      0.195    
3       Aruba vs Denmark  1  4.54372380 50.8948430  0.637022880   0.001      0.015   .
4      Aruba vs Malaysia  1  0.87672195  6.0397240  0.172367910   0.013      0.195    
5         Aruba vs Spain  1  3.49955224 36.9984243  0.560595570   0.001      0.015   .
6    Brazil vs CaboVerde  1  1.79573937 27.7696018  0.489163230   0.001      0.015   .
7      Brazil vs Denmark  1  0.27026868 35.0545409  0.547260825   0.001      0.015   .
8     Brazil vs Malaysia  1  1.69391038 26.6369446  0.478763612   0.001      0.015   .
9        Brazil vs Spain  1 -0.01258568 -0.9666596 -0.034482498   0.953      1.000    
10  CaboVerde vs Denmark  1  2.85656348 47.6207158  0.629731090   0.001      0.015   .
11 CaboVerde vs Malaysia  1 -0.02164970 -0.1836835 -0.006603446   0.978      1.000    
12    CaboVerde vs Spain  1  1.56212224 23.8545897  0.460028512   0.001      0.015   .
13   Denmark vs Malaysia  1  2.82610262 48.0023359  0.631590271   0.001      0.015   .
14      Denmark vs Spain  1  0.38029019 58.5480612  0.676480332   0.001      0.015   .
15     Malaysia vs Spain  1  1.47580731 22.9256154  0.450178466   0.001      0.015   .

#pdf("figures/beta_q1f_loca.pdf",width=9, height=5)
beta_q1f$S %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = "sample") %>%
  group_by(location) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = location, fill = location)) +
    geom_point(size = 4) +
    #   stat_ellipse(aes(color = beta_q1n_nmds$Groups))+
    geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
    scale_color_manual(values = location_colors)+
    theme_classic() +
    theme(
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      axis.title = element_text(size = 20, face = "bold"),
      axis.text = element_text(face = "bold", size = 18),
      panel.background = element_blank(),
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      legend.text = element_text(size = 16),
      legend.title = element_text(size = 18),
      legend.position = "right", legend.box = "vertical"
    ) 

#dev.off()

6.2 Behaviour

6.2.1 Richness diversity

betadisper(beta_q0n$S, sample_metadata$origin) %>% permutest(., pairwise=TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.00411 0.0041147 0.2442    999  0.613
Residuals 90 1.51634 0.0168483                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Domestic Feral
Domestic          0.621
Feral     0.62238      

adonis2(beta_q0n$S ~ origin, 
        data = sample_metadata %>% arrange(match(sample,labels(beta_q0n$S))), 
        permutations = 999,
        strata = sample_metadata %>% arrange(match(sample,labels(beta_q1n$S))) %>% pull(location)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
origin	1	0.4407951	0.0163738	1.498173	0.178
Residual	90	26.4799587	0.9836262	NA	NA
Total	91	26.9207538	1.0000000	NA	NA

pairwise.adonis(beta_q0n$S, sample_metadata$origin, perm = 999)

              pairs Df SumsOfSqs  F.Model        R2 p.value p.adjusted sig
1 Domestic vs Feral  1 0.4407951 1.498173 0.0163738   0.081      0.081    

beta_q0n$S %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(origin) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = origin, fill = origin)) +
    geom_point(size = 4) +
    scale_color_manual(values = origin_colors) +
    geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
    theme_classic() +
    theme(
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      axis.title = element_text(size = 20, face = "bold"),
      axis.text = element_text(face = "bold", size = 18),
      panel.background = element_blank(),
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      legend.text = element_text(size = 16),
      legend.title = element_text(size = 18),
      legend.position = "right", legend.box = "vertical"
    ) 

6.2.2 Neutral diversity

betadisper(beta_q1n$S, sample_metadata$origin) %>% permutest(., pairwise=TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.0080 0.0080048 0.5581    999  0.459
Residuals 90 1.2909 0.0143435                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Domestic Feral
Domestic          0.463
Feral     0.45698      

adonis2(beta_q1n$S ~ origin, 
        data = sample_metadata %>% arrange(match(sample,labels(beta_q1n$S))), 
        permutations = 999,
        strata = sample_metadata %>% arrange(match(sample,labels(beta_q1n$S))) %>% pull(location)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
origin	1	0.386668	0.01498735	1.369385	0.315
Residual	90	25.412962	0.98501265	NA	NA
Total	91	25.799630	1.00000000	NA	NA

pairwise.adonis(beta_q1n$S, sample_metadata$origin, perm = 999)

              pairs Df SumsOfSqs  F.Model         R2 p.value p.adjusted sig
1 Domestic vs Feral  1  0.386668 1.369385 0.01498735    0.14       0.14    

pdf("figures/beta_q1n_behaviour.pdf",width=9, height=5)
beta_q1n$S %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(origin) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = origin, fill = origin)) +
    geom_point(size = 4) +
    scale_color_manual(values = origin_colors) +
    geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
    theme_classic() +
    theme(
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      axis.title = element_text(size = 20, face = "bold"),
      axis.text = element_text(face = "bold", size = 18),
      panel.background = element_blank(),
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      legend.text = element_text(size = 16),
      legend.title = element_text(size = 18),
      legend.position = "right", legend.box = "vertical"
    ) 
dev.off()

RStudioGD 
        2 

6.2.3 Phylogenetic diversity

betadisper(beta_q1p$S, sample_metadata$origin) %>% permutest(., pairwise=TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.01229 0.012295 0.5442    999  0.464
Residuals 90 2.03337 0.022593                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Domestic Feral
Domestic          0.467
Feral     0.46262      

adonis2(beta_q1p$S ~ origin, 
        data = sample_metadata %>% arrange(match(sample,labels(beta_q1p$S))), 
        permutations = 999,
        strata = sample_metadata %>% arrange(match(sample,labels(beta_q1n$S))) %>% pull(location)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
origin	1	0.2191073	0.02185403	2.010807	0.16
Residual	90	9.8068363	0.97814597	NA	NA
Total	91	10.0259435	1.00000000	NA	NA

pairwise.adonis(beta_q1p$S, sample_metadata$origin, perm = 999)

              pairs Df SumsOfSqs  F.Model         R2 p.value p.adjusted sig
1 Domestic vs Feral  1 0.2191073 2.010807 0.02185403   0.073      0.073    

beta_q1p$S %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(origin) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = origin, fill = origin)) +
    geom_point(size = 4) +
    scale_color_manual(values = origin_colors) +
    geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
    theme_classic() +
    theme(
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      axis.title = element_text(size = 20, face = "bold"),
      axis.text = element_text(face = "bold", size = 18),
      panel.background = element_blank(),
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      legend.text = element_text(size = 16),
      legend.title = element_text(size = 18),
      legend.position = "right", legend.box = "vertical"
    ) 

6.2.4 Functional diversity

betadisper(beta_q1f$S, sample_metadata$origin) %>% permutest(., pairwise=TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.2486 0.248604 3.3502    999  0.088 .
Residuals 90 6.6786 0.074206                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Domestic Feral
Domestic          0.091
Feral    0.070508      

adonis2(beta_q1f$S ~ origin, 
        data = sample_metadata %>% arrange(match(sample,labels(beta_q1f$S))), 
        permutations = 999,
        strata = sample_metadata %>% arrange(match(sample,labels(beta_q1n$S))) %>% pull(location)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
origin	1	0.3816603	0.02530502	2.336579	0.18
Residual	90	14.7007312	0.97469498	NA	NA
Total	91	15.0823914	1.00000000	NA	NA

pairwise.adonis(beta_q1f$S, sample_metadata$origin, perm = 999)

              pairs Df SumsOfSqs  F.Model         R2 p.value p.adjusted sig
1 Domestic vs Feral  1 0.3816603 2.336579 0.02530502   0.127      0.127    

beta_q1f$S %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(origin) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = origin, fill = origin)) +
    geom_point(size = 4) +
    scale_color_manual(values = origin_colors) +
    geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
    theme_classic() +
    theme(
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      axis.title = element_text(size = 20, face = "bold"),
      axis.text = element_text(face = "bold", size = 18),
      panel.background = element_blank(),
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      legend.text = element_text(size = 16),
      legend.title = element_text(size = 18),
      legend.position = "right", legend.box = "vertical"
    )